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.1  Introduction 

The  advantage  of  the  far  ultraviolet  radiation  in  the  investigation  of  airglow  and 
aurora  is  that  pure  absorption  is  strong  enough  to  prevent  the  radiation  from 
penetrating  to  depths  where  rayleigh  scattering,  ground  and  cloud  albedo,  and 
scattering  from  aerosols  can  take  place. 

The  USAF  Geophysics  Laboratory  has  developed  and  flown  the  first  instru¬ 
mentation  designed  for  UV  imagery  with  high  spatial  and  wavelength  resolution. 
The  AIM  scanner  flown  on  the  HILAT  satellite  provided  the  earliest  data  show¬ 
ing  airglow  and  auroral  morphology  and  emission  processes.  HILAT  returned 
a  number  of  images  in  1983  with  a  30A  spectral  resolution  and  5  x  20  km  spa¬ 
tial  resolution.  The  AIM  instrumentation  aboard  HILAT  has  been  described  in 
several  references  ([5],  [8],  [11],  and  [17]). 

In  late  1986,  AFGL  launched  a  second  satellite  with  similar  but  improved 
instrumentation  ([2],  and  [18]).  The  Polar  BEAR  satellite  has  been  designed 
for  UV  imagery  with  high  spatial  and  wavelength  resolution.  The  satellite  car¬ 
ries  the  AIRS  scanning  ultraviolet  photometer  providing  several  simultaneous 
spectral  images.  These  images  represent  the  integrated  column  photon  count 
of  several  different  upper-atmosphere  species.  The  data  stream  is  currently 
converted  to  images  representing  5000  X  5000  Jkm  records  of  intensity. 

The  records  of  intensity  obtained  by  the  Polar  BEAR  satellite  become  mean¬ 
ingful  results  only  by  manipulating  them  by  image  processing  techniques. 

Previous  research  by  our  group,  originally  led  by  the  late  Prof.  N.  Rosenberg, 
has  laid  the  foundations  for  the  photometric  and  geometric  processing  of  both 
the  AIM  and  the  AIRS  instruments  ([12],  [13],  [14],  [15],  and  [  16})- 

A  software  package  has  been  developed  to  process  large  number  of  images 
supplied  by  the  AIRS  instrument.  The  package  transforms  the  time-spaced 
data  stream  to  a  photometrically  and  geometrically  standardized  image  in  a 
semiautomatic  manner. 

Table  1  presents  a  summary  of  programs  and  procedures  newly  developed 
or  updated  in  the  coarse  of  this  research.  A  data-flow  diagram  for  basic  image 
processing  of  the  AIRS  pictures  is  shown  in  Figure  1.  The  flowsheet  shows 
processes,  and  the  relations  between  them.  A  single  scene  is  defined  as  a  multi¬ 
detector  record  for  a  given  site,  day,  and  pass. 
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name 

description 

amov 

automatic  image  moving  procedure 

(screen  to  screen,  screen  to  file,  and  vise  versa) 

attitude 

fit  an  offset  sinsuidal  function  8 (t)  to  attitude  angles; 
extrapolate  6  to  fill-in  null  pitch,  yaw,  and  roll 

aurlaq 

search  for  auroral  brightness  correlation  with  path  length 

aunnean 

average  multiple  auroral  images;  produces  multiple 
image-coverage,  auroral  frequency,  and  average  intensity 

aurplot 

plot  image  upon  corrected  geomagnetic  coordinate  system 

chapman 

Chapman  function  Ch(x,  X)  for  65°  <  *  <  105°,  X  =  170 

circle 

draws  a  circle  of  chosen  color  and  radius  on  screen 

circum 

plot  image  circumference  (field-of-view) 

compraw 

compare  two  orbital  data  files 

corrgm 

corrected  geomagnetic  lat,  Ion  for  geographic  lat,  Ion 

craur 

fit  a  circle  to  maximum  auroral  brightness 

detc 

mode-smooth,  amplify,  filter,  and  warp  image 

dipgm 

dipole  geomagnetic  lat,  Ion  for  geographic  lat,  Ion 

distance 

distance  image  of  object  to  background  and  vise  versa 

edgim 

find  smallest  rectangle  bounding  image 

frame 

sketch  border  lines  on  screen 

geomlsq 

correlate  dayglow  intensity  with  path  length  and  x; 
subtract  fitted  background  dayglow  from  warped  image 

Table  1:  Basic  Software  Package 
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name 


description 


gm  coord 

overlay  geographic  or  geomagnetic  grid  upon  image 

grafaur 

graph  auroral  intensity  versus  magnetic 
latitude,  local-time,  or  pixel  distance 

grafile 

graph  any  ascii  dth  integer  pairs  on  screen 

grafim 

graph  image  histogram,  or  (average)  line/column  brightness 

ht 

halftone  image  from  byte  to  bit 

ht2hp 

printer  utility  -  halftone  to  laser  printer  (hp) 

ht2pr 

printer  utility  -  halftone  to  matrix  printer  (epson) 

imfair 

get  Interplanetary  Magnetic  Field  data  for  AIRS  image 

integrate 

integrate  image  brightness 

landmans 

overlay  geographic  grid  and/or  landmass  features  on  image 

latoval 

plot  polar  corrected  geomagnetic  latitude — magnetic 
local-time  coordinate  system  and/or  Feldstein’s  oval 

Iinlsq 

linear  least-squares  fit  to  n  data  points  for  m  unknowns 

list 

list  orbital  data  summary  of  multiple  orbital  files 

med 

calculate  attitude  parameters  statistics 

medfill 

automatic  procedure  for  filling-in  image  with  median  filter 

monthday 

month  and  day  from  year  and  day  number 

mulbnd 

two-color  output  image  from  two  input  detectors 

newscl 

re-scales  image  brightness 

overiim 

overlay  one  image  on  another 

Table  1:  Basic  Software  Package  (cont) 
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name 

description 

plattit 

graph  attitude  parameter  pitch,  yaw,  or  roll  versus  time 

pllat 

automatic  procedure  to  plot  geomagnetic  latitude  lines 

plotgm 

draw  dipole  or  corrected  geomagnetic  grid  upon  image 

plovlat 

interactively  plot  coordinate  system  and  auroral  oval 

pLgraf 

automatic  procedure  to  graph  two  integer  vectors  on  screen 

rawcvt 

convert  raw  ascii  file  to  binary  image  +  orbital  data  file 

replsq 

report  least-squares  fit  data  of  multiple  images 

rewarp 

interactive  program  to  re-warp  raw  image 

rhline.p 

emulate  PC  integer  size  to  fiV ax  integer  size 

rhline.v 

emulate  FG100AT  procedures  to  IP-512  driver  software 

scene 

calculate  orbital,  geometrical,  and  geophysical  parameters 

stdsube 

standard  subroutines,  including  Euler’s  rotation  matrix 

sundec 

i 

sun  declination  from  date  and  time 

terminator 

mark  sun  terminator  on  an  image  mapped  by 
either  geomagnetic  or  geographic  coordinates 

warpm 

automatic  procedure  to  warp  an  image  with  error  control 

xmov 

interactive  image  moving  program 

(screen  to  screen,  screen  to  file,  and  vise  versa) 

zoomfll 

zoom  image  file  to  image  file  (in  multiples  of  2) 

loomim 

magnify  image  on  screen 

Table  .1:  Basic  software  package,  (cont) 


4 


Figure  .1:  Data-flow  diagram  of  AIRS  image  processing. 
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.2  Pre-processing 

.2.1  Raw  Data  Structure 

The  Polar  BEAR  satellite  cruises  at  an  altitude  of  approximately  1000  kilome¬ 
ters  at  an  inclination  angle  of  89.56  degrees  to  the  equator,  i.e.  in  a  nearly  polar 
orbit. 

The  Auroral  Ionospheric  Remote  Sensor  (AIRS)  is  the  prime  instrument  on¬ 
board  the  three-axis  stabilised  craft.  AIRS  provides  326  samples  of  emission 
count  in  an  angular  scan  from  -65.2°  to  +65.2°  perpendicular  to  the  satel¬ 
lite  orbit.  The  scan  cycle  is  completed  in  3  seconds  and  a  typical  11  minutes 
recording  provides  ~  220  scan  lines. 

AIRS  acquires  a  single  scene  per  orbit  at  the  two  working  wavelengths 
(1356.4,  and  1596A).  Each  pixel,  which  is  equivalent  to  a  mirror  step,  has 
an  integration  time  of  6.83  milliseconds,  and  a  readout  time  of  0.2  msec. 

The  VUV  Instantaneous  Field  Of  View  (IFOV)  is  0.373°  x  1.53°.  Raw  image 
structure  is  summarised  in  Figure  2,  in  which  each  ellipse  is  a  single  pixel.  Line 
and  column  are  l  and  Jfc  respectively. 


Figure  .2:  Raw  VUV  image  structure. 
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The  VUV  pixels  are  mapped  on  an  emitting  layer  150  km  above  ground. 


There,  the  nadir  pixels  have  the  following  dimensions: 

a 

= 

0.373° 

=  5.5  km 

b 

= 

1.53° 

=  22.7  km 

c 

= 

0.4° 

=  5.9  km 

d 

= 

1.3° 

=  19.3  km 

At 

= 

t|+i  - 1| 

--  3  sec 

Ar 

= 

Tk+l  -  Tk 

=  7  msec 

The  spacecraft  forward  velocity  o,  was  designed  to  be  6.6  km  /sec  at  ground 
subtract  point.  We  fitted  ground  path  d  to  time  of  flight  t  via 

d  =  d0  +  v,  ■  (t  -  t0) 

for  each  scene  to  find  v,  =  6.3  ±0.02  km/sec  for  all  scenes  processed  so  far.  This 
velocity  brings  about  a  ground  distance  of  18.9  km  between  consecutive  lines, 
and  a  distance  of  19.3  km  between  consecutive  lines  at  an  altitude  of  150  km. 

Because  the  satellite  platform  is  not  well  stabilized,  pitch,  yaw,  and  roll 
corrections  to  the  scan  pointing  direction  must  be  included.  In  general  scan 
lines  are  not  parallel,  and  are  not  perpendicular  to  the  orbital  ground  path. 

.2.2  Binary  Image  Preparation 

Each  raw  scan  record  contains  two  lines  —  a  data  line  of  326  valuable  intensity 
values,  and  a  control  line  with  the  following  data: 

•  Universal  Time  (UT)  of  scan  initiation  (in  seconds). 

•  Geographic  latitude  and  longitude  of  satellite  nadir. 

•  Satellite  attitude  angles:  pitch,  yaw,  and  roll. 

•  Height  of  satellite  above  ground  (in  kilometers). 

The  initial  processing  (program  RAWCVT)  converts  a  raw  ascii  file  with 
records  as  above  to  two  files: 

1.  A  binary  image  of  240  rows  x  240  columns. 

2.  An  ascii  file  containing  41  lines  of  orbital  data. 
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The  orbital  data  file  contain*  control  lines  from  the  original  file  for  every 
sixth  record.  The  image  is  formed  by  converting  ~  220  lines  of  numbers  to  240 
lines  of  binary  data.  The  number  of  lines  in  a  raw  file  varies  from  120  to  330, 
with  an  average  of  220.  Lines  are  either  multiplied  or  skipped  upon  in  order  to 
conform  the  binary  image  to  the  predetermined  site.  The  warping  of  the  image 
in  the  post-processing  stage  eliminates  the  redundancy  in  identical  lines,  as  will 
be  shown  later.  The  case  of  line  skipping  is  extremely  rare,  as  most  raw  files 
contain  240  records  or  less.  The  program  handles  data  gaps  and  tests  whether 
the  raw  data  is  packed.  If  packed,  procedure  UNPAK  is  invoked  to  generate  a 
nonsaturated  image. 

The  sise  of  the  binary  image,  a  matrix  of  240  X  240  pixels,  is  maintained  in 
order  to  conform  to  past  work  and  data  files.  The  resampling  of  the  columns  is 
more  complicated  than  line  stretching.  This  part  of  the  processing  has  gone  a 
major  change  since  past  work,  which  influences  the  image  structure  considerably. 
Hence  it  will  be  discussed  separately  in  the  following  subsection. 

.2.3  Resampling  the  Raw  Line 

The  raw  line  contains  336  values,  of  which  the  first  and  last  five  are  null.  These 
two  unsampled  degrees  on  each  side  of  the  scan  line  correspond  to  the  fact  that 
the  mirror  scans  through  ±67.2  degrees,  of  which  only  ±65.2  degrees  of  data  is 
collected  due  to  telemetry  constraints. 

Figure  3  shows  the  resampling  of  326  pixels  in  a  raw  line  to  240  pixels  of  the 
output  line: 


output  line 


Figure  .3:  Resampling  geometry. 

The  first  and  last  meaningful  14  pixels  (column  #  5-18,  and  317-330)  of  the 
raw  line  are  copied  directly  to  the  first  and  last  14  pixels  (column  #  0-13,  and 
226-239)  of  the  output  line  respectively.  These  28  pixels  retain  their  signal  and 
IFOV  perpendicular  to  flight  path — 0.4°.  The  inner  298  pixels  of  the  input  line 
are  resampled  to  the  inner  212  pixels  of  the  output  line. 

The  differential  sampling  is  due  to  the  fact  that  we  want  to  maintain  maxi¬ 
mum  resolution  at  the  edges  of  the  image.  We  use  these  edge  pixels  later  on  to 
correct  the  roll  of  the  satellite  by  adjusting  it  to  the  location  of  the  limb. 
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At  the  central  part  of  the  raw  line,  we  transform  each  input  pixel  with  an 
IFOV  of  0.4°,  to  an  output  pixel  with  an  IFOV  of  ■  0.4°  ~  0.56°. 

The  poor  signal  in  each  pixel  drives  us  to  maintain  the  raw  line-integrated 
intensity.  Therefore  we  redistribute  the  signal  of  the  inner  298  pixels  to  212 
pixels  with  a  normalised  weighting  function. 

The  redistribution  scheme  is  according  to  weights  proportional  to  the  dis¬ 
tance  of  the  designated  index  to  the  two  closest  integers.  Let  g  be  the  dis¬ 
crete  weighting  function,  *  the  input  index,  and  j  the  output  index.  Then  for 
*  =  19,. ..,316: 

x  =  (t  —  19)  •  +  14 

ii  =  L*J 

ft  =  1-  (*-  l*J) 


If  IN  ia  the  input  line,  and  OUT  the  output  line,  then 

316 

OUTu  =  'E'W  .ft 

♦=i» 

316 

OUTh+l  =  £  IN,  (1  -  e.) 

•  =  19 

A  major  deficiency  of  the  above  weighted  redistribution  is  its  nonlinearity, 
i.e.  the  sum  of  weights  per  output  pixel  is  not  exactly  constant.  The  fact  that  the 
transfer  function  is  nonlinear  implies  that  the  resampling  of  a  constant  signal , 
i.e.  a  line  with  a  constant  value,  is  not  constant. 

A  linearisation  method  is  introduced  in  order  to  normalise  the  sum  of  weights 
of  every  output  index  to  the  actual  increase  in  IFOV  (jjj|  =  1.406).  If  g  is  the 
desired  normalised  weighting  function,  then 

316 

fa  =  Y, « 

•=10 

316 

fai+i  =  EU-ft) 

|  =  19 

and  for  j  =  14, ... ,  225: 

fa  =  1-406/  fa 

Once  the  output  line  has  been  computed  based  on  the  input  line,  the  intensity 
of  each  output  pixel  is  normalised  to  the  actual  increase  in  IFOV,  via 

OUTj  =  OUTj  fa 
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Figure  4  shows  a  typical  OUT  image,  as  produced  by  RAWCVT. 


Note  the  relatively  low  signal  that  is  present  in  the  unstretched  image.  Count 
dynamic  range  is  50,  equivalent  to  a  maximum  intensity  of  ~  IK  Rayleighs 
(1  connt  =  21  Rayleighs  for  a  pixel  with  an  IFOV  of  0.56°  at  1356A). 


The  raw  image  is  customarily  rotated  to  align  north  up,  east  to  the  right 
side  of  the  image.  Thin  auroral  arcs  inhabit  the  central  part  of  the  image. 
The  scattered  radiation  at  the  lower  half  of  the  image  is  twilight  and  dayglow. 
Nightglow  intensity  is  below  the  threshold  of  the  detector,  hence  no  counts  are 
registered  in  the  upper  mid-part  of  the  image.  Note  the  brightened  limbs,  to  be 
discussed  in  the  following  section. 


Figure  .4:  Raw  image  -  file  s87023h2  (scene  07),  detector  2  (A  =  1 356 A) . 
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.2.4  Attitude  Control 

The  Polar  BEAR  is  a  three-axis  stabilised  satellite.  A  son  sensor  and  a  flux  gate 
magnetometer  measure  the  three  angles  of  the  instantaneous  position  of  the 
satellite  relative  to  the  nadir  direction.  Since  one  of  the  craft  attitude  sensors 
is  a  sun  sensor,  attitude  data  is  available  only  under  sunlight. 

On  the  average,  the  instantaneous  instability  of  the  spacecraft  is  of  order  1-2 
degrees  in  pitch,  2-5  degrees  in  roll,  and  4-11  degrees  in  yaw.  These  magnitudes 
vary  with  time  on  a  time-scale  of  weeks  to  months.  Over  the  entire  lifetime  of 
the  satellite,  these  instabilities  have  caused  the  satellite  to  invert  to  the  senith 
direction  more  than  once. 

At  night,  no  attitude  data  is  available.  Almost  all  images  are  partly  in 
the  earth’s  shadow.  Thus,  in  principle,  each  and  every  image  can  be  partially 
correctly  oriented,  partially  not.  This  fact  creates  a  distortion  not  only  between 
images,  but  also  within  a  single  image. 

The  null  hypothesis  is  that  the  satellite  points  directly  downward  at  the  nadir 
direction.  This  allows  one  to  assume  sero  attitude  corrections  to  the  three  di¬ 
rections  in  non-sunlit  conditions.  In  light  of  the  above  magnitudes  of  pitch,  roll, 
and  especially  yaw,  as  measured  in  the  sunlit  hemisphere,  we  find  this  limita¬ 
tion  quite  significant  when  one  tries  to  find  the  geographic  and  geomagnetic 
coordinates  of  the  image  correctly. 

In  order  to  be  able  to  correct  for  the  satellite  attitude  variations  when  no 
corresponding  data  is  available,  we  have  embarked  on  the  following: 

•  Error  quantification  using  a  few  typical  images. 

•  Study  of  the  Polar  BEAR’s  attitude  changes  in  time. 

•  Development  of  an  algorithm  to  solve  the  problem,  in  effect  stabilising  the 
craft  retrospectively. 

•  Test  the  statistical  validity  of  the  technique. 

Error  quantification 

Obviously,  the  error  in  position  is  minimal  in  the  nadir  direction,  and  maximal 
in  the  east  and  west  edges  of  each  line. 

We  scaled  two  images,  obtained  under  full  sunlight,  with  and  without  the 
attitude  corrections  as  provided  by  the  sun  sensor  and  the  magnetometer.  The 
images  were  chosen  to  represent  data  with  close  to  minimum  (t8703072 — scene 
24),  and  close  to  maximum  (r87051k2 — scene  40)  attitude  corrections.  Pitch, 
roll,  and  yaw  amplitudes  in  scene  24  were  2°,  3°,  and  4s  respectively,  and  2°, 
4°,  and  10°  respectively  in  scene  40. 
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Table  2  shows  ground  distances  x  and  y  (in  kilometers)  at  the  edges  of  the 
two  pairs  of  data.  The  cartesian  x,  y  coordinate  system  is  laid  on  a  hypothetical 
spherical  earth.  The  x-axis  is  the  projection  of  the  satellite  track  on  that  sphere; 
the  y-axis  is  perpendicular  to  it.  The  y  coordinate  is  formed  by  the  projection  of 
the  instantaneous  line-of-sight  to  the  center  of  the  emitting  layer  (assumed  to  be 
at  an  altitude  of  150  km) ,  and  down  to  the  earth.  The  origin  of  the  coordinate 
system,  i.e.  (x,  y)  =  (0,0),  is  at  the  nadir  location  of  the  satellite  at  first-line 
acquisition  time. 


without  attitude  corrections 

scene  24 

with  attitude  corrections 

1  *  ..  -  . 

y _ 

*  _ 

y 

2 

0 

-4 

-1968 

0 

1960 

-49 

-15 

5 

MM 

-19 

1968 

1852 

1849 

1844 

-1982 

0 

1982 

1942 

1824 

1680 

Hi 

-14 

1976 

3702 

3699 

3695 

•  1998 

0 

1998 

3698 

3685 

3659 

-15 

1997 

scene  40 

without  attitude  corrections 

with  attitude  corrections 

x 

y 

x  ..  . 

_ y  1 

1 

-1 

-5 

-2014 

0 

2014 

-408 

-25 

335 

-1978 

1 

1990 

2104 

2101 

2096 

-2023 

0 

2023 

2002 

2070 

2108 

-2020 

13 

2023 

4206 

4203 

4199 

-2028 

0 

2028 

4390 

4175 

3935 

-2021 

0 

2013 

Table  .2:  Influence  of  attitude  corrections  on  ground  distances  (in  km). 

For  the  best  case,  i.e.  scene  24,  the  averaged  central  x  and  y  distances  differ 
by  18  and  16  kilometers,  corresponding  to  ~  1,  and  ~  3  VUV  pixels  respectively. 
At  the  edges,  x  and  y  distances  without  attitude  corrections  differ  from  the 
distances  that  include  attitude  corrections  by  as  much  as  47  and  70  kilometers 
respectively. 

By  contrast,  scene  40  exhibits  differences  of  up  to  200  kilometers  in  the 
x-axis  at  the  edges,  and  similar  errors  in  the  central  pixels.  The  greater  error 
in  the  x  distances  relative  to  the  y  distances  is  due  to  a  large  yaw  correction . 

The  conclusion  is  therefore  that  the  error  in  pixel  localisation  when  using  an 
image  lacking  attitude  data  is  significant  even  for  cases  where  pitch,  roll,  and 
yaw  are  2-5°,  and  even  at  the  central  part  of  the  image.  If  large  (>  10°)  attitude 
corrections  are  in  order,  the  error  per  pixel  is  of  order  100-300  kilometers  at  the 
edges,  and  tapers  off  to  a  few  dosens  of  kilometers  at  the  center. 

Temporal  attitude  behavior 

Rosenberg  and  Quesada  ([16])  noted  that  the  Polar  BEAR  satellite  changes  its 
attitude  in  a  periodic  fashion.  We  have  found  this  periodicity  to  sustain  (with 
increasing  amplitude  in  time)  over  the  first  half  year  of  Polar  BEAR’s  life  in 
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orbit  by  plotting  a  few  dosen  passes  randomly  chosen  from  that  era. 

Figure  5  shows  the  temporal  evolution  of  the  three  attitude  angles  for  scene 
40,  as  extracted  from  file  rS7051k2: 


/ 
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Figure  .5:  Attitude  angles  versus  time  -  scene  40. 


The  pitch,  yaw,  and  roll  angles  of  the  satellite  vary  periodically  with  time. 
The  temporal  behavior  of  the  three  attitude  angles  can  be  approximated  by  an 
offset  sine  function: 


9(f)  =  Ae  sin(iue  t  +  <f>e)  +  ©(0) 

where  the  attitude  angle  0,  as  a  function  of  time  t,  depends  on  four  variables: 
Ae  -  amplitude,  -  frequency,  <f>s  -  phase,  and  8(0)  -  offset. 

By  fitting  all  four  parameters  independently  to  several  scenes,  we  found  that 
the  function  describes  the  three  angles  well  Hence  it  seems  reliable  to  describe 
the  motion  of  the  satellite  around  three  orthogonal  axes  by  above  function. 


Periodic  attitude  extrapolation 

The  amplitude,  phase,  frequency,  and  offset  magnitudes  change  slowly  in  time 
(on  a  time-scale  of  weeks  to  months).  Hence  we  devised  the  following  algorithm 
to  extrapolate  daytime  attitude  corrections  to  nighttime.  The  method  can  be 
applied  within  a  certain  scene,  or  between  two  different  scenes  closely  spaced  in 
time. 

An  offset  sinusoidal  is  fitted  to  attitude  data  as  reported  by  the  two  sensors 
in  the  day  side  part  of  the  image.  The  extrapolated  coefficients  are  used  to  fill-in 
the  missing  attitude  data  for  the  nighttime  part  of  the  scene.  The  extrapolation 
for  a  scene  with  very  little  or  no  attitude  data  at  all  is  performed  based  on 
data  from  another  scene  with  sufficient  attitude  information.  Sufficient  data  is 
defined  to  be  a  data  set  with  a  time  span  larger  than  half  a  period  (P@  =  ^). 
The  scene  from  which  the  parameters  are  established  from  should  be  temporally 
close  to  the  designated  scene  as  possible,  and  never  more  than  a  few  days  apart. 

Figure  6  shows  an  application  of  the  self-extrapolation  technique  on  the  yaw 
data  (program  ATTITUDE).  Bar  histogram  line  shows  original  data,  smooth 
line  is  the  fitted  and  extrapolated  one.  The  plot  itself  was  produced  by  PLATTIT. 
Scene  54  (file  s8705992)  was  used  to  produce  the  results  shown  in  Figure  6  and 
the  following  table.  In  this  least-squares  fit  a  free  choice  of  frequency  for  every 
attitude  angle  was  allowed.  The  resultant  periods  are  395.2,  390.3,  and  395.2 
seconds  for  yaw,  roll,  and  pitch  respectively.  These  results  support  Rosenberg 
and  Quesada’s  (|16j)  statement  regarding  the  coupling  between  the  three  angles. 
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Figure  .6:  Yaw  data  extrapolation  -  scene  54. 


attitude 

angle 

-^e 

amplitude 

(degrees) 

Pe 

period 
( seconds ) 

phase 

(degrees) 

0(0) 

offset 

(degrees) 

yaw 

7.5 

395.2 

51.8 

0.4 

roll 

1.0 

390.3 

-53.1 

-2.3 

pitch 

0.6 

395.2 

-123.8 

1.2 

Table  .3:  Determined  variables  of  the  attitude  self-extrapolation  technique. 


Statistical  tests 

We  took  the  following  steps  in  order  to  test  the  alidity  of  the  above  technique: 

•  Choose  two  scenes,  closely  spaced  in  time,  both  taken  under  sunlit  condi¬ 
tions;  the  two  scenes  then  contain  complete,  or  almost  complete  attitude 
information. 

•  Use  one  scene’s  data  to  extrapolate  pitch,  roll,  and  yaw  to  the  second  scene 
and  vise-versa;  the  result  will  be  two  pairs  of  data  sets,  corresponding 
to  the  two  scenes;  each  pair  includes  the  actual  attitude  angles  and  the 
extrapolated  ones. 
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•  Perform  a  statistical  comparison  of  each  pair  to  find  the  correlation  coef¬ 
ficient  between  the  actual  and  the  calculated  attitude  angles. 

Ideally  the  two  sets  of  data  should  match,  and  we  should  get  a  perfect 
correlation  (  r  =  1  )  between  the  two  sets.  If  the  extrapolated  values  differ  from 
the  actual  ones  significantly,  the  coefficient  of  correlation  r  should  be  low  or 
insignificant  according  to  the  t  test. 

Table  4  presents  the  results  of  a  test  of  this  nature  performed  by  program 
COMPRAW  between  scene  06  and  scene  07.  The  dates  and  times  of  the  two 
scenes  are:  scene  06  •  year  87  -  day  47  -  start  time  43229.15  seconds,  and  scene 
07  -  year  87  -  day  48  -  start  time  51323.15  seconds. 

File  1  contains  the  raw  data,  file  2 — the  calculated  data.  Notation  of  statis¬ 
tical  parameters  is  as  follows:  fie  -  arithmetic  mean,  9e  -  mean  deviation,  cr© 
-  standard  deviation,  100- -  coefficient  of  variation,  5e,|e,  -  standard  error 
of  the  estimate,  r  -  coefficient  of  correlation,  and  t  -  t  test. 


scene  0  file#  fie  oe  100  ^  Se,ie 


06  roll 


07  roll 


1  0.76  5.03  5.69  113.07% 

2  1.06  5.37  6.03  112.26% 

1—2 


1  -3.00  0.67  0.77  25.73 

2  -2.53  0.74  0.82  32.33% 

1—2 


1  1.17  0.27  0.33 

2  1.99  0.28  0.31  15.71% 

1—2 


1  1.27  5.19  5.91  113.96% 

2  1.51  4.75  5.43  114.34% 

1—2 

1  -2.64  0.77  085  32.09% 

2  -2.87  0.64  0.70  24.48% 

1—2 


1  1.97  0.44  0.52  26.16% 

2  1.14  0.24  0.27  23.54% 


1.53  0.97  23.73 


0.57  0.92  14.70 


0.85  0.77  7.27 


2.56  0.90  12.03 


0.41  0.92  13.67 


0.94  0.50  3.33 


Table  .4:  Actual  and  extrapolated  attitude  angles  -  a  statistical  comparison. 

Note  that  t  >  ;.5  implies  a  level  of  significance  P  of  null  hypothesis  (0i  and 
©2  are  uncorrelated)  of  0.001,  hence  the  value  of  r  is  significant! 

The  results  of  the  above  test  clearly  indicate  an  excellent  agreement  between 
the  actual  parameters  and  the  extrapolated  ones.  The  value  of  the  coefficient  of 
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correlation  is  above  90%  for  yaw  and  roll.  It  is  lower  for  the  pitch  angle  since 
its  amplitude  is  relatively  smalL  It  is  comparable,  in  fact,  to  the  1-2°  error 
measurement  of  these  angles  by  the  attitude  sensors. 


.3  Geometric  Processing 

.3.1  Motivation 

Geometric  corrections  are  applied  to  the  image  in  order  to  rectify  and  map  it 
onto  the  ionosphere.  Furthermore,  satellite  attitude  corrections  are  imperative 
to  assure  that  image  coordinates  are  properly  rotated  to  the  real-time  body 
reference  system.  To  the  best  knowledge  of  the  authors,  this  project  is  the  first 
one  to  fully  include  geometric  and  photometric  corrections  to  VUV  images  of 
the  ionosphere. 

The  geometric  corrections  include  the  following: 

•  A  correction  to  the  satellite  roll  based  on  the  dayglow  limbs  observed  at 
the  edges  of  the  image. 

•  Image  rectification  via  a  transformation  which  is  a  function  of  the  instan¬ 
taneous  satellite  attitude  angles. 

•  The  geographic  location  of  each  pixel  within  the  image  is  located  by  tracing 
the  satellite  line-of-sight  to  the  height  of  the  emitting  layer. 

.3.2  Roll  Adjustment 

Evident  from  the  raw  image  shown  in  Figure  4,  the  east  and  west  edges  of  the 
image  are  typified  by  an  intensity  enhancement. 

Assume  that  the  ionospheric  region  contributing  to  the  airglow  intensity 
has  a  constant  thickness.  If  the  emission  is  optically  thin,  then  the  above 
phenomenon  is  easily  explained  by  limb  brightening.  Even  if  some  multiple 
scattering  effects  cause  this  phenomenon,  it  is  still  clear  that  it  is  path  length 
dependent. 

Hence  we  receive  an  increased  radiation  as  the  mirror  scans  from  local  nadir 
to  the  horison.  All  other  parameters  constant,  the  maximum  intensity  will  be 
at  the  point  where  the  length  of  the  emitting  layer,  relative  to  the  detector,  is 
maximum. 

Since  not  all  the  other  parameters  equal,  the  limb  enhancement  occurs  un¬ 
evenly  for  the  two  sides  of  the  image.  In  Figure  4,  for  example,  the  sun  is  located 
toward  the  southeast  edge  of  the  image,  making  the  east  (right)  limb  brighter 
than  the  west  (left)  one. 

Based  on  the  above,  we  use  the  brightest  limb  that  is  present  in  the  image 
to  correct  for  the  satellite  roll.  Roll  adjustment  is  performed  in  three  steps: 
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1.  The  brightest  limb  in  the  image  is  chosen  according  to  its  angular  prox¬ 
imity  to  the  sun. 

2.  Search  for  the  brightest  column  in  every  image  line  in  the  neighborhood 
of  tne  chosen  edge.  A  linear  interpolation  is  performed  in  order  to  £11- 
in  column  numbers  for  unsuccessful  lines.  The  search  itself  is  under  a 
continuity  restraint,  i.e.  column  numbers  in  successive  lines  can  differ  by 
2  pixels  at  the  most. 

3.  The  adjustment  to  satellite  roll  per  :mage  line  is  the  difference  between 
the  scan  angle  that  corresponds  to  the  column  number  of  maximum  in¬ 
tensity,  and  the  scan  angle  that  corresponds  to  a  line-of-sight  tangent  to 
the  emitting  layer. 

The  roll  adjustment  is  performed  only  on  limbs  under  dayglow  radiation, 
i.e.  for  solar  senith  angles  x  <  88°.  The  theoretical  scan  angle  of  line-of-sight 
tangent  to  the  center  of  the  emitting  layer  is  geometrically  computed  using  a 
Chapman  approximation  to  limb  height  h: 

h  =  160  -  25  ■  log10(cos  \) 

see  also  Strickland  and  Anderson  (|20|). 

The  following  table  lists  the  average  roll  adjustment  of  a  few  scenes  with 
corresponding  solar  senith  angles  (x)  at  the  center  of  the  image. 


scene 

# 

X  at 

image  center 

roll 

adjustment 

06 

88.0 

1.9  ±  0.3 

07 

92.4 

1.5  ±  0.2 

32 

96.0 

1.8  ±  0.2 

33 

96.4 

1.9  ±  0.3 

34 

80.5 

1.5  ±  0.3 

35 

84.3 

1.7  ±  0.4 

36 

86.5 

1.6  ±  0.8 

37 

84.6 

1.5  ±  0.3 

38 

87.7 

1.4  ±  0.4 

40 

90.5 

1.2  ±  0.4 

57 

89.2 

1.2  ±  0.4 

59 

90.6 

to 

H- 

o 

Table  .5:  Average  roll  adjustment. 


For  a  specific  scene,  the  roll  may  vary  by  as  much  as  3°,  whereas  the  roll 
adjustment  has  a  dynamic  range  of  ~  1°. 

Since  something  like  a  minimum  intensity  criteria  (  x  <  88°)  is  imposed  on 
the  brightness  of  the  limb  in  order  to  initiate  this  procedure,  it  is  clear  that 
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images  obtained  at  local  midnight  will  have  no  roll  adjustment.  Major  parts  of 
the  north  polar  region  experience  continuous  darkness  during  the  winter.  That 
implies  that  quite  a  few  images  at  that  time  will  lack  the  roll  correction. 

Notice,  however,  that  the  average  roil  adjustment  per  scene,  as  listed  in 
table  5,  is  quite  constant.  The  variation  of  the  average  roll  adjustment  between 
scenes  is,  in  fact,  much  smaller  than  both  the  variation  of  the  roll  and  the  roll 
adjustment  within  a  certain  scene.  It  may  be  that  the  attitude  sensors  have  a 
constant  offset  of  its  readings.  Thus  we  apply  a  constant  roll  adjustment  of  1.6° 
in  case  the  roll  adjustment  procedure  may  not  be  applied. 

.3.3  Basic  Geometry 

Figure  7  shows  satellite  -  earth  geometry  as  used  in  this  project.  The  satellite 
cruises  at  an  altitude  of  H,  kilometers  above  ground.  The  ionospheric  emitting 
layer  responsible  for  the  intensity  measured  by  the  satellite  lies  between  altitudes 
/ij  and  h<2}  with  thickness  Ah  =  hi  —  h 2.  We  use  h  =  <  hj ,  h2  >  as  the  actual 
height  of  the  emitting  layer  for  pixel  localisation. 

Satellite  reference  system  is  represented  by  the  (X,,  Y„  Z.)  coordinate  sys¬ 
tem.  X,  is  the  roll  axis,  Y,  the  pitch  axis,  and  Z,  -  the  yaw  axis.  Scan  angle  is 
measured  from  Z,  on  the  Y,Z,  plane. 

The  reference  coordinate  system  is  (X,Y,Z).  Its  X-axis  is  directed  toward 
v,,  the  satellite  velocity  vector.  The  Z-axis  originates  within  the  satellite  and  is 
directed  toward  earth’s  center.  The  y-axis  complements  a  right-handed  system. 
The  reference  coordinate  system  follows  the  motion  of  the  satellite  in  its  orbit. 

In  general,  one  wishes  to  find  the  coordinates  of  pixels  within  the  image  in 
a  geocentric  coordinate  system.  Two  transformations  are  required  : 

•  From  the  body  (satellite)  coordinate  system  to  a  reference  system. 

•  FYom  the  reference  system  to  a  geocentric  system,  which  is  (almost)  inertial. 


Vector  l,  in  the  ( X,,Y,,Z ,)  coordinate 
system,  points  toward  the  instantaneous 
line-of-sight.  It  has  an  angle  a  to  Z,,  a 
being  the  scan  angle.  Vector  p,  in  the 
(X,  F,  Z)  coordinate  system,  has  an  angle 
of  a'  degrees  to  the  Z  axis.  It  is  a 
transformation  of  /  via  : 

p  =  A  f 


Figure  .7:  Satellite  -  Earth  geometry. 
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where 


P  = 


co*  a' 


2  =  1  sin  a 
\  coe  a 


Following  Sommerville  ([19]),  A  =  (a,,)  is  Euler’s  general  rotation 

matrix  : 


an 

= 

R2 

cosO  +  (1-cosfi) 

“12 

= 

Y  PR 

q  sin fl  +  n3  (1  cosO) 

“13 

= 

--  smfl  +  (l-cosfl) 

“21 

= 

^  sinft  +  (1  cosO) 

“22 

= 

P2 

cosO  +  —  (1-cosfl) 

«»23 

= 

R  .  _  P  Y 
-  smO  +  -jp-  ■  (1  —  cos  fl) 

“31 

= 

^  sinfl  +  ~  (1-cosO) 

“32 

= 

n  sinO  +  (1  cosn) 

“33 

= 

Y2 

cos  n  +  ( i  -  cos  n) 

with  n2  =  R2  +  P2  +  Y2  ,  and  R,  P,  and  Y  as  the  roll,  pitch,  and  yaw 
angles  respectively. 

The  length  of  the  vector  p,  i.e.  p,  is  the  distance  between  the  satellite  and  the 
emitting  layer  at  an  altitude  of  h  kilometers  above  ground  for  the  instantaneous 
scan  angle  a  : 


p  =  (Rq  +  H „)  cos  o'  -  yj (R9  +  H ,)2  cos2  a'  -  ( R e  +  H,)2  +  (R9  4-  h)2 

The  vertex  of  the  angle  ft  is  O,  earth’s  center.  It  is  the  angle  between  -Z 
and  the  vector  from  O  to  the  point  at  which  p  is  tangent  to  the  emitting  layer  : 


f)  =  arcsin 


(  v !pl+pl\ 
^  R9  +  h  J 
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.3.4  Image  Ground  Coordinates 

The  (XeiYe)  system  (see  Fig.  7)  is  a  2-D  cartesian  coordinate  system  on  the 
surface  of  a  spherical  earth.  The  Xe  axis  is  the  projection  of  the  satellite  track 
on  the  ground.  The  Ye  is  perpendicular  to  it.  The  origin  of  the  system  is  at  the 
ground  projection  of  the  satellite  location  at  first-line  acquisition  time. 

The  coordinates  Xe  and  Ye  are  ground  distances.  Yc  is  the  distance  of  the 
projection  of  the  scan  line-of-sight  to  satellite  nadir.  Xe  is  the  ground  distance 
of  the  instantaneous  satellite  nadir  point  to  the  position  of  the  first  line  parallel 
to  the  projected  satellite  track. 

Yr  and  Xt  are  given  by  the  solution  of  a  spherical  triangle  : 


X,.  =  R$  ■  arctan(cos  T  •  tan  0) 

Yt  =  if®  ■  arcsin(sin  T  ■  sin/?) 


where  T  is  the  angle  between  the  projection  of  p  on  the  X  —  Y  plane  and  the 
X  axis. 

The  geographic  latitude  (<j>)  and  longitude  (A)  of  the  projection  of  each  pixel 
on  the  emitting  layer  can  be  computed  via  its  Xe,  Ye  coordinates.  Let 


0  = 


n. 

R& 


and  get  <f>  and  A  as  a  function  of  0  and  cr: 


sin  t  •  sin  <7  +  0  cost 

Vl  +  03 

tan  Ao  ■  cos  a  +  cos  *  •  sin  a  —  ©  •  sin  « 
cos  a  —  tan  Ao  ■  cos  i  ■  sin  a  +  0  •  tan  Ao  •  sin  » 

The  inclination  of  the  orbit  is  i.  Ao  is  the  geographic  longitude  of  the  orbit 
at  the  equator.  <ro  is  the  angular  distance  (on  the  ground)  between  the  orbital 
intersection  with  the  equator  and  the  origin  of  Xc,Ye. 


sin  <t>  = 

tan  A  = 


.3.5  Image  Rectification 

The  ground  mapping  of  pixels  on  the  Xe,Ye  coordinate  system  is  performed 
only  for  certain  scan  angles.  A  scan  angle  must  be  in  the  range  of  the  east  and 
west  scan  angles  that  are  tangent  to  the  spherical  ground  earth.  Other  scan 
angles  are  ambiguous  in  its  ground  locations,  since  its  line-of-sights  crosses  the 
center  of  the  emitting  layer  more  than  once. 

The  mapping  is  performed  on  41  equally  spaced  scan  angles.  Thus  a  41  X  41 
matrix  is  mapped  on  the  surface  of  the  earth.  Each  4  point  cell  of  the  above 
matrix  is  used  for  a  bi-linear  interpolation  of  all  36  points  of  the  original  image 
within  that  cell. 
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Since  the  attitude  of  the  satellite  affects  the  scan  direction,  the  scan  in 
general  in  neither  perpendicular  to  the  orbit,  nor  symmetric  relative  to  the 
nadir.  Moreover,  the  variable  yaw  causes  the  scan  lines  to  be  non-parallel. 

Hence  the  image  must  be  rectified.  We  desire  an  image  with  lines  perpen¬ 
dicular  to  the  satellite  orbit  and  columns  parallel  to  it.  A  transformation  is 
created  from  the  raw  image  to  the  warped  image.  Each  pixel  of  the  raw  image 
is  located  by  its  ground  Xe,Ye  coordinates  using  the  technique  above.  Then 
its  line  and  column  number  in  the  rectified  image  are  calculated  via  the  closest 
neighbor  method. 

The  method  of  closest  neighbor  employs  the  distance  of  a  warped  pixel  to 
the  closest  grid  point.  The  intensity  of  a  rectified  grid  point  is  the  intensity  of 
the  original  pixel  with  minimum  distance  to  that  point.  The  distance  d  is  given 
by 

d  =  n/A*2  +  AP 


where  A k  and  A 1  are  the  differences  between  warped  column  and  line  and 
nearest  grid  column  and  line  respectively. 

Figure  8  shows  a  distance  image  produced  by  WARPM  procedure.  Increas¬ 
ing  distance  is  marked  by  an  increased  blackening  of  the  image. 


Figure  .8:  Distance  image  -  scene  59  (file  r87052l2). 

If  the  entire  set  of  raw  pixels  is  to  fall  within  the  borders  of  the  warped 
image,  an  appropriate  scale  must  be  found  for  the  new  map.  The  chosen  sise  of 
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the  warped  image  is  4800  km  x  4800  km.  Consequently,  some  edge  pixels  within 
the  rectified  image  have  no  inverse  transformation  to  the  raw  image,  and  are 
thus  null.  The  empty  border  regions  are  marked  by  the  grey  tone  in  Figure  8. 

The  transformation  from  the  raw  image  to  the  rectified  image  gives  us  an 
opportunity  to  change  the  way  the  spherical  earth  is  projected  onto  the  two 
dimensional  image.  The  imaging  method  employed  by  the  satellite  creates  a 
space  oblique  mercator  projection  of  the  ionospheric  north  pole  region.  The 
transformation  of  the  raw  pixels  using  a  fixed  absolute  distance  ( =  20  ^7^7) 
creates  a  polar  stereographic  projection. 

.4  Photometric  Processing 

.4.1  Motivation 

The  photometric  processing  enhances  the  image  in  order  to  make  use  of  the 
full  dynamic  range  afforded  by  the  image  processing  hardware.  In  a  subsequent 
development,  we  devised  a  novel  approach  toward  the  separation  of  dayglow 
from  aurora.  The  technique  is  based  on  a  spatially  variable  intensity  criteria. 
The  dayglow  intensity  in  every  single  image  is  fitted  by  a  general  function. 
Variables  of  this  function  are  the  solar  senith  angle,  and  the  path  length  of  the 
emitting  layer.  Pixels  with  an  intensity  larger  than  that  function  are  defined  as 
auroral  ones. 


.4.2  Geophysical  Parameters 

Four  geophysical  variables  are  computed: 

•  Inverse  path  length. 

•  Corrected  cos  x- 

•  Number  of  hours  since  sunrise  or  sunset. 

•  Total  daily  solar  exposure  time. 

These  4  geophysical  parameters  are  stored  in  a  single  binary  image  240  x  240 
pixels  large.  It  may  be  displayed  as  an  ordinary  image,  where  each  parameter 
occupies  an  area  of  120  lines  X  120  columns  X  1  byte.  Figure  9  shows  this  image 
for  scene  06  (file  s87023f2).  Inverse  path  length  is  on  the  upper-right,  cosx  on 
the  upper-left,  daily  solar  time  on  the  lower-right,  and  time  since  sunrise/sunset 
on  the  lower-left  side  of  the  image. 

The  path  length  D  of  the  emitting  layer  as  observed  by  the  detector  for  a 
specific  scan  angle  a  is  given  by  D  =  Pht  —  PH,,  where  pnt  and  p*,  are  the 
lengths  of  the  vector  p  as  it  crosses  the  emitting  later  at  heights  h  1  and  hj 
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respectively.  Since  A h  <  H®,  D  can  be  approximated  by 


D  _  _ (1?$+  h) _ 

AA  y/(Rm  +  H.f  cos3  <*'-(*•  +  if.)2  +  {R*  +  h )a 

The  solar  senith  angle  x  depends  on  the  declination  of  the  sun  S,  the  hour 
angle  t  of  the  sun,  and  the  geographic  latitude  ^  : 

cos  x  =  sin  6  sin  <f>  +  cos  S  ■  cos  <j>  ■  cos  t 

The  sun’s  declination  angle  is  given  by  ainS  —  sin  e  •  sin  l,  where  e  ~  23.5°  is 
the  inclination  of  the  ecliptic  to  the  equator,  and  l  the  celestial  longitude  of  the 
sun,  i.e.  the  longitude  of  the  sun  on  the  ecliptic  relative  to  the  equinox. 

The  amount  of  solar  radiation  hitting  the  atmosphere  affects  the  airglow 
intensity.  However,  the  functional  dependence  of  cos  x  is  incorrect  at  high  solar 
senith  angles  due  to  the  general  failure  of  the  plane-parallel  approximation  at 
these  angles  ((lj).  By  high  solar  senith  angles  we  mean  90°+  if>  >  x  >  65°, 
where  if>  is  the  150  km  UV  sunset/sunrise  dip  angle.  We  correct  cos  x  using 
modified  Chapman  function  CA(x,  X)  values  with  a  scale  height  gradient  of 
=  0  75  at  150  km.  X  =  =  170,  where  H{0)  =  38.4  km  is  the 

atomic  oxygen  scale  height  at  150  Am  ([21],  and  [23]).  Following  Kasachevskaya 
et  al.  ([10]),  we  assume  if>,unTi,e  =  10°,  and  V'.un.rt  =  5°. 


Figure  .9:  Four  geophysical  parameters  •  selOO.i 


Solar  radiation  function 


Figure  10  shows  the  difference  between  Ch(x,X)  and  cos  x  as  a  function 


Figure  .  10:  A  comparison  of  Ch(x,  X )  and  cos  x  f°r  100°  >  ^  >  65°. 
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.4.3  Enhancement  and  Filtering 

As  evident  from  Figure  4,  the  dynamic  range  of  the  raw  image  is  rather  poor. 
Detector  2  (A  =  1596A)  has  a  dynamic  range  of  ~  50,  whereas  the  dynamic 
range  of  detector  1  (A  =  1356A)  is  ~  5.  In  order  to  utilise  the  full  dynamic 
range  afforded  by  the  image  processing  hardware  (255),  a  stretching  routine  is 
applied  to  the  image.  The  amplification  factor  m  is  based  on  the  particular 
image  statistics.  It  is  set  to  m  =  255/6gg.g,  where  igg.g  is  the  brightness  of  the 
99.9  percentile  in  the  cumulative  histogram  of  the  raw  image. 

Since  the  signal  to  noise  ratio  of  the  image  is  limited  by  count  statistic 
itignai  _  l  \  the  low  signal  implies  a  relatively  high  noise  level.  The  effect 

'  *****  yj ttgnal 

of  non-discriminate  linear  amplification  on  noisy  pixels  is  highly  unwelcome. 
Hence  a  smoothing  routine  is  applied  to  the  amplified  image.  We  use  the  most- 
frequent  value  filter  in  order  to  eliminate  single  noise  pixels. 

If  the  image  features  some  airglow  background,  it  may  be  desired  to  subtract 
it.  The  subject  of  background  subtraction  will  be  discussed  in  the  following 
subsection.  However,  the  background  airglow  must  be  smoothed  before  trying 
to  fit  any  continuous  function  to  it.  Thus  we  apply  two  averaging  filters  to  the 
image  before  initiating  the  background  subtraction:  an  auroral  ,  and  an  airglow 
averaging  filter.  The  auroral  filter  has  the  dimensions  of  5  columns  x  3  lines; 
the  airglow  filter  is  15  columns  x  7  lines  large.  We  have  chosen  the  ratio  of  filter 
width  by  height  to  be  ~  2  to  approximately  compensate  for  the  angular  ratio 
of  the  distance  between  resampled  pixels  (1.3°/0.56°). 

.4.4  Background  Subtraction 

The  night  glow  intensity  is  below  the  detection  threshold  of  the  two  working 
detectors  (  1,  and  2).  Therefore  when  one  discusses  airglow  subtraction  in  Polar 
BEAR  images,  one  actually  refers  to  dayglow  and  twilight  glow  only. 

Airglow  in  general  is  defined  in  a  negative  sense.  It  is  the  entire  atmospheric 
radiation  excluding  aurora,  Rayleigh  scattering,  and  thermal  emission.  The  sole 
source  of  airglow  energy  is  solar  radiation.  Dayglow  consists  of  prompt  radiation, 
and  night  glow  of  delayed  emission.  Twilight  glow  is  dayglow  observed  from  the 
night  sky.  Thus  the  transmission  of  twilight  glow  differs  from  the  transmission 
of  dayglow  through  the  atmosphere. 

However,  not  enough  twilight  glow  pixels  are  present  in  a  typical  image  to 
warrant  a  statistically  valid  fit.  Hence  a  single  function  is  fitted  to  the  two 
airglow  domains.  Our  aim  here  is  to  fit  an  analytic  function  of  geophysical 
parameters  to  the  geometrically  corrected  image  in  order  to  separate  aurora 
from  airglow. 

A  least-squares  fit  is  performed  on  the  brightness  b  of  a  pixel  as  a  function 
of  two  variables:  path  length  l,  and  the  Chapman  function  s  (corrected  cos(x))- 
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The  fit  of  the  function 


b  =  c0  •  pc’  ■  «f  * 

produces  three  coefficients:  coi  cpi  and  c,. 

Results  of  such  fits  to  9  scenes  and  2  detectors  are  presented  in  Table  6.  c0 
is  the  average  brightness  (eo)  normalised  to  unit  path  length  (i  =  1)  and  sero 
senith  distance  (s  =  1).  uo  is  the  standard  deviation  around  the  mean  of  log10  6, 
where  b  is  the  original  pixel  brightness,  oj  is  the  standard  deviation  around  the 
mean  of  the  expression  log10  b  —  cp  log10/  —  c,log10  a.  Ratios  ^  <  1  indicate 
the  successfulness  of  the  fit. 


scene  j 
# 

M 

mmm 

HI 

HI 

06 

100.4 

0.70 

1.08 

0.18 

0.11 

89.8 

0.77 

0.84 

0.16 

0.03 

07 

111.6 

0.69 

1.02 

0.15 

0.08 

221.9 

0.67 

1.19 

0.13 

0.04 

34 

81.0 

0.70 

0.98 

0.19 

0.10 

77.5 

0.82 

1.03 

0.19 

0.04 

35 

45.7 

0.76 

0.60 

0.19 

0.11 

64.3 

0.83 

0.88 

0.18 

0.05 

37 

87.5 

0.60 

0.90 

0.18 

0.11 

63.9 

0.84 

0.82 

0.19 

0.06 

38 

68.5 

0.58 

0.43 

0.17 

0.12 

66.8 

0.84 

0.75 

0.19 

0.07 

40 

26.0 

1.10 

0.46 

0.23 

0.08 

62.6 

0.88 

0.87 

0.20 

0.03 

57 

39.6 

0.90 

0.47 

0.20 

0.10 

57.4 

0.94 

0.79 

0.22 

0.04 

59 

45.8 

0.62 

0.70 

0.14 

0.06 

103.2 

0.69 

1.12 

0.18 

0.04 

Table  .6:  Least-squares  fit  to  airglow  intensities. 

Program  GEOMLSQ  uses  the  coefficients  found  to  subtract  the  expected 
from  the  actual  brightness  of  every  pixel.  Following  this  subtraction,  if  the 
brightness  of  a  pixel  is  larger  than  the  amplification  factor  m,  then  it  is  an 
auroral  one.  A  subsidiary  criteria  of  an  auroral  brightness  threshold  is  used 
in  conjunction  with  the  above  subtraction  in  order  to  distinguish  airglow  from 
auroral  pixels. 

The  program  uses  the  warped  image  and  the  geophysical  parameters  to 
create  an  auroral  image.  An  image  that  compares  the  fitted  intensity  values 
with  the  actual  ones  is  also  provided. 

Once  a  purely  auroral  image  has  been  formed,  the  raw  image  is  rewarped  in 
order  to  produce  a  non- averaged  auroral  image.  Program  REWARP  re-warps 
the  raw  image  while  checking  the  auroral  image  as  produced  by  GEOMLSQ. 
Non-auroral  pixels  are  set  to  a  signal  of  1  (pixels  in  the  warped  image  that  carry 
a  sero  brightness  are  pixels  outside  the  imaging  field-of-view). 

Figure  11  shows  the  same  image  shown  in  Figure  4  photometrically  and 
geometrically  processed.  Most  of  the  noise  has  been  smoothed  out,  and  the 
majority  of  airglow  background  subtracted. 
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Figure  .11:  Auroral  image  -  scene  07,  A  =  1356A. 


.5  Geometric  Error  Analysis 

The  major  purpose  of  the  geometric  processing  is  to  map  the  image  correctly 
on  the  geographic  coordinate  system.  The  location  of  each  pixel  is  a  function 
of  the  satellite  attitude  parameters  and  the  height  of  the  emitting  layer.  The 
accuracy  of  the  mapping  depends  on  the  accuracy  of  these  parameters. 

We  have  analysed  the  effect  of  the  main  uncertainties,  i.e.  satellite  attitude 
angles  and  height  of  emitting  layer,  on  the  geographic  location  of  pixels  as  a 
function  of  scan  angle. 

The  accuracy  of  the  control  system  is  specified  as  1°.  However,  our  roll  cor¬ 
rection  technique  implies  an  error  of  1-2°  for  the  roll  angle,  and  by  implication 
for  the  pitch  and  yaw  angles  as  well.  We  have  assumed  a  constant  height  of 
150  km  for  the  emitting  layer,  and  believe  its  accuracy  to  be  ~  50  km. 

A  statistically  random  approach  has  been  chosen  for  the  error  analysis.  A 
large  number  of  iterations  were  run  to  compute  the  geographic  location  of  pixels 
for  normally  random  distributions  of  the  above  variables.  The  result  was  a  dis¬ 
tribution  of  the  mapped  location,  for  which  the  average  and  standard  deviations 
were  computed. 

The  standard  deviation  from  the  exact  location  (sero  errors  in  variables)  is 
given  by  an  angle  whose  apex  is  at  the  center  of  the  earth.  One  may  assume, 
then,  that  the  maximum  error  in  either  latitude  or  longitude  (geographic  or 
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geomagnetic)  is  this  angle. 

Let  P  be  the  pitch  angle,  R  the  roll  angle,  and  Y  -  the  yaw  angle.  <  P>, 
<  R  >,  and  <  Y  >  are  the  average  values,  and  op,  or,  and  ay  the  standard 
deviations  of  these  angles  respectively.  We  used  these  six  values  plus  the  height 
of  the  emitting  layer  ( H  =  150  ±  50i:m)  to  create  random  groups  of  these 
variables  normally  distributed  around  the  mean.  We  chose  <P>  =  1.7°,  <R> 
=  —1.0°,  and  <Y>  =  7°.  The  standard  deviations  are  op  =  ay  =  1-2°,  with 
the  exception  of  or  =  1°  due  to  the  roll  correction. 

The  simulation  was  run  for  seven  scan  angles:  0° ,  30°,  40°,  50°,  55°,  60° ,  and 
62°  for  three  cases  : 

•  Error  in  attitude  angles  only. 

•  Error  in  emitting  layer  height  only. 

•  Error  in  both  attitude  angles  and  emitting  layer  height. 

Table  7  shows  the  results  of  these  error  simulation  runs.  Each  case  was 
iterated  5000  times  to  create  a  the  statistic  of  a  normal  distribution.  To  test 
the  randomness  of  the  procedure,  the  average  values  and  standard  deviations  of 
the  distributions  of  attitude  angles  were  compared  to  the  corresponding  input 
values.  6  is  the  average  value  of  the  error;  at  -  the  standard  deviation  of  the 
error  (in  degrees). 
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Table  .7:  Angular  error  in  pixel  ground  location. 
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The  error  increases  with  scan  angle.  It  is  larger  than  1°  for  scan  angles  >60°. 
In  the  worst  case,  an  average  error  of  2°  exists  for  a  pixel  located  at  a  scan  angle 
of  60°. 

The  results  of  the  error  simulations  also  point  to  a  practical  aspect.  We  can 
give  an  upper  limit  of  1°  on  the  error  anywhere  in  the  image  if  we  limit  the 
mapping  to  scan  angles  relative  to  satellite  nadir  that  are  less  than  55  degrees. 
An  error  of  1°  is  quite  acceptable,  since  it  is  also  the  accuracy  of  the  corrected 
geomagnetic  coordinate  system  ([9j),  and  the  accuracy  of  Feldstein’s  oval  ([22|), 
to  be  discussed  in  the  following  section. 

Thus  a  calculation  of  a  lower-limit  scan  angle  derived  from  an  angle  of  55° 
relative  to  the  nadir  is  performed  by  program  SCENE.  In  the  course  of  pixel 
projection  to  its  ground  location,  procedure  WARPM  utilises  this  limit  to  avoid 
mapping  pixels  beyond  the  corresponding  scan  angle.  The  price  paid  is  loss  of 
a  small  part  of  the  original  image  coverage,  and  converting  linear  image  edges 
to  edges  that  follow  the  satellite  roll. 

We  compare  the  mapping  of  an  image  without  and  with  the  above  procedure 
in  the  following  two  Figures.  Figure  12  shows  an  auroral  map  in  which  all 
possible  scan  angles  are  used.  Figure  13  depicts  the  same  scene  (scene  44  •  file 
s87053l2)  where  the  mapping  has  been  limited  to  scan  angles  that  correspond 
to  angles  <  55°  relative  to  the  nadir. 


Figure  .12:  Mapping  of  an  image  without  a  55°  nadir  angle  limit. 
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Figure  .13:  Mapping  of  an  image  with  a  55°  nadir  angle  limit. 

.6  Coordinate  Systems 

.6.1  Geographic  System 

Two  programs  are  available  for  overlaying  a  geographic  grid  upon  an  image: 
LANDMASS  and  GMCOORD.  LANDMASS  plots  a  geographic  grid  and/or 
landmass  features  upon  the  warped  image.  GMCOORD  plots  a  geographic  grid 
upon  the  image  with  a  somewhat  better  accuracy  than  LANDMASS. 

A  note  should  be  made  regarding  the  landmass  feature  of  the  first  program. 
The  pixels  of  the  warped  image  are  mapped  onto  a  sphere  150  kilometers  above 
ground,  whereas  landmasses  are  on  the  ground  earth.  The  image  and  Iandmasses 
are  thus  projected  to  a  common  picture  from  a  view  point  of  satellite  nadir  at 
image  central  line. 

Program  TERMINATOR  marks  the  sun  terminator  on  a  240x240  image 
mapped  by  cartesian  geographic  coordinates.  A  user  may  choose  the  jolar  dip 
angle  of  the  terminator.  The  default  value  is  set  to  the  150  kilometers  sunset 
solar  dip  angle  (-12.3°). 

Figure  14  shows  scene  06  (file  s87023f2)  for  which  no  airglow  has  been  sub¬ 
tracted.  A  geographic  grid  has  been  overlaid  by  GMCOORD  and  landmass 
features  were  drawn  by  LANDMASS.  Two  terminators  are  marked:  the  geo¬ 
metric  horison  at  solar  senith  angle  \  =  90°,  and  the  150  km  sunset  /sunrise 
terminator  at  x  =  102.3°. 
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Figure  .14:  Image  062  with  a  geographic  grid,  landmasses,  and  terminators. 


.6.2  Geomagnetic  System 

Two  different  techniques  are  available  for  geomagnetic  coordinate  system  appli¬ 
cation.  GMCOORD,  as  in  the  geographic  case,  plots  a  geomagnetic  grid  upon 
the  image.  AURPLOT  does  exactly  the  opposite  -  it  warps  the  image  onto  a 
fixed  coordinate  system.  Both  allow  the  use  of  two  different  coordinate  systems: 
magnetic  dipole,  and  corrected  geomagnetic. 

lhe  magnetic  dipole  is  a  first  approximation  to  the  geomagnetic  field.  The 
imaginary  dipole  is  located  at  the  center  of  the  earth,  and  inclined  by  11.5° 
relative  to  earth’s  rotation  axis.  The  dipole  crosses  the  earth’s  surface  at  78.5° N, 
291.0° E,  which  is  the  north  pole  of  the  geomagnetic  dipole  system. 

The  field  lines  are  given  by  r  =  ro  •  sin2  0;  ro  is  the  geocentric  distance 
where  the  field  line  crosses  the  magnetic  equator. 

Note  that  the  dipole  system  is  simply  a  spherical  system  similar  to  the  geo¬ 
graphic  system.  Thus  any  point  on  the  earth  with  geographic  coordinates  4>,  A 
is  translated  to  the  geomagnetic  dipole  system  by  rotating  the  system  around 
the  geographic  axis  by  291  degrees  and  around  the  new  axis  by  11.5  degrees. 


The  coordinate  system  most  used  in  the  literature  is  the  corrected  geomag¬ 
netic  system.  Hakura  ([7j)  computed  the  ground  shift  of  the  magnetic  field 
lines  due  to  the  first  four  terms  in  the  spherical  harmonic  development  of  the 
potential  following  the  dipole  term.  Hakura’s  system  was  improved  by  Gustafs- 
son  ((6j),  who  included  two  more  harmonic  terms.  Gustafsson  showed  that  the 
maximum  contribution  of  these  terms  is  130  kilometers,  to  which  the  7-th  term 
contributed  less  than  30  km. 

The  corrected  system  is  computed  by  following  a  common  point  on  the 
equator  along  a  dipole  field  line,  and  along  a  corrected  field  line  (with  the 
higher  terms)  to  the  ground.  The  ground  distance  on  the  sphere  is  then  the 
correction  to  the  dipole  system  for  that  specific  point. 

Gustafsson  created  a  table  to  convert  geographic  coordinates  to  corrected 
geomagnetic  coordinates  by  using  the  first  seven  harmonic  terms.  The  table 
supplies  corrected  geomagnetic  coordinates  for  every  2°  of  geographic  latitude, 
and  every  5°  of  geographic  longitude.  Intermediate  points  are  calculated  by 
interpolation. 

The  corrected  system  of  Gustafsson  is  a  ground  system,  whereas  the  aurora 
occurs  at  an  altitude  of  —150  kilometers.  A  dipole  correction  to  the  corrected 
geomagnetic  latitude  9  c  as  a  function  of  height  h  is  given  by 

cos2  9 c(h)  =  R<*o+  h  cos2  9C(0) 

A$ 

where  R9  is  earth’s  radius. 

The  warping  of  an  image  to  a  fixed  geomagnetic  polar  stereographic  system 
leaves  ’holes’,  similar  to  the  ones  created  when  magnifying  an  image.  These 
holes  are  filled  by  a  median  filter. 

The  approach  of  projecting  an  image  to  a  fixed  coordinate  system  affords 
us  to  view  and  compare  all  the  images  on  a  common  ground.  Furthermore, 
multiple-image  operations  are  readily  available. 

Figure  15  shows  two  images  taken  during  a  single  pass  over  the  north  pole 
by  two  different  receiving  stations.  Image  402  (file  r87051k2)  and  412  (file 
t87051k2)  are  rotated  to  align  local  magnetic  noon  at  the  top  of  the  map. 
Corrected  magnetic  latitude  lines  are  marked  at  10°  intervals;  magnetic  local 
time  at  6  hours  intervals. 

The  two  images  span  a  total  integration  time  of  22  minutes,  of  which  2 
minutes  are  overlapping  time.  The  102.3°  terminator  of  this  overlapping  time  is 
the  arc  at  the  central  part  of  the  map.  Background  airglow  has  been  subtracted. 
The  sigsag  border  is  the  overall  field-of-view  of  the  two  images.  The  limit  to 
scan  angle  of  55°  relative  to  satellite  nadir  has  been  applied. 
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Figure  .15:  Images  402  and  412  on  a  corrected  magnetic  time-latitude  system. 
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.7  Results  and  Conclusions 

.7.1  Multiple  Auroral  Image  Averaging 

Feldstein  et  al.  ([3],  and  [4])  prepared  ovale  of  auroral  activity  as  a  function 
of  the  geomagnetic  Q  index.  Feldstein’s  ovals  are  the  results  of  optical  images 
obtained  by  all  sky  cameras  on  the  ground.  Within  the  oval  boundaries  there 
is  a  chance  of  70%  for  the  aurora  to  occur  at  any  time.  The  boundaries  are  the 
results  of  averaging  the  poleward  and  equatorward  boundaries  of  the  aurora  as 
detected  by  these  cameras.  The  ovals  represent  mainly  discrete  aurora,  since 
the  continuous  aurora  exhibits  a  much  lower  contrast,  and  is  much  weaker  in 
intensity  per  square  sky  degree  than  the  discrete  aurora. 

The  fixed  coordinate  system  described  in  the  previous  section  allows  us  to 
average  many  images,  and  create  the  VUV  version  of  Feldstein’s  ovals.  A  major 
advantage  this  work  has  on  Feldstein’s  is  the  fact  that  we  can  include  daylight 
aurora.  Also  the  relatively  small  bandwidth  of  the  detector  spectral  sensitivity 
(~  30 A)  implies  we  are  looking  mainly  at  a  single  species  (OI). 

Program  AURMEAN  combines  auroral  images  to  produce  three  maps  -  im¬ 
ages  coverage,  normalised  auroral  frequency,  and  average  intensity.  Normalised 
frequency  is  the  number  of  times  aurora  was  present  in  a  certain  piece  of  sky 
relative  to  the  total  number  of  times  that  piece  of  sky  was  imaged. 

We  have  processed  45  scenes  at  two  wavelengths:  1596A  (detector  1),  and 
1356A  (detector  2).  The  images  were  grouped  according  to  integer  values  of 
the  Kp  index,  and  the  Q  index.  The  binning  was  required  in  order  to  gather 
a  reasonable  number  of  images  per  group  for  averaging.  Also  Feldstein’s  ovals 
are  grouped  according  to  Q  values,  so  Kp  binning  is  a  must  for  this  kind  of 
comparison.  Images  with  Kp  =  2+,  3,  and  3-,  for  example,  were  binned  under 
the  Kp  =  3  group. 
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Table  .8:  Scene  numbers  according  to  Kp  and  Q  values. 


Table  8  lists  all  scenes  that  contributed  to  the  VUV  ovals.  A  single  scene 
inhabits  the  Kp  =  5  group,  therefore  ignored.  Although  the  sise  of  the  entire 
sample  is  reasonable,  considering  the  amount  of  processing  required,  the  sample 
per  group  is  rather  small  This  effects  both  the  degree  to  which  the  value  is 
representative  of  the  value  of  a  large  distribution  of  averaged  images,  and  the 
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total  coverage  of  the  images. 

Figures  16-25  depict  auroral  frequency  as  a  function  of  Kp.  The  local  mag¬ 
netic  time  -  corrected  geomagnetic  latitude  system  is  used.  Borders  describe  the 
total  field  of  view  of  the  ensemble  of  images.  Notice  that  in  almost  all  cases  the 
total  coverage  of  the  averaged  images  does  not  supply  us  with  a  full  coverage  of 
the  upper  part  of  the  northern  geomagnetic  hemisphere. 

Figures  26-35  depict  average  auroral  intensity  as  a  function  of  Kp.  The  same 
coordinate  system  is  used.  The  ovals  of  Feldstein  are  present  in  Figures  16-35. 

The  conclusions  to  be  drawn  from  the  VUV  ovals  are  as  follows  : 

•  The  aurora  in  ultraviolet  radiation  appears  in  ovals  that  surround  the 
geomagnetic  pole;  the  center  of  the  ovals  is  not  at  the  pole. 

•  The  latitude  span  of  the  aurora  is  greatest  around  local  geomagnetic  mid¬ 
night;  it  is  4°  at  Kp  =  4. 

•  As  the  magnetic  activity  rises  { Kp  increases),  the  poleward  boundary  of 
the  aurora  turns  from  an  oval  to  a  circle  with  a  radius  of  ~  14°  centered 
on  the  geomagnetic  pole. 

•  In  general  there  is  a  good  correlation  between  the  VUV  ovals  and  Feld- 
stein’s  oval;  a  conspicuous  exception  is  the  dayside  aurora  at  low  Kp;  at 
Kp  =  0,  1  the  maximum  latitude  difference  is  5°. 

There  are  a  few  plausible  explanations  to  the  discrepancy  between  our  maps 
and  Feldstein 's  ovals  at  low  Kp  around  local  magnetic  noon.  Since  the  statistical 
coverage  of  the  dayside  is  rather  poor  both  in  our  study  and  in  Feldstein ’s,  it  may 
simply  be  a  statistical  error.  We  intend  to  expand  our  processing  in  the  future 
to  include  dosens  of  images  in  each  map.  Thus  this  point,  at  least  regarding 
the  statistical  validity  of  our  maps,  will  most  certainly  be  cleared  in  the  near 
future. 

One  major  difference  that  exists  between  Feldstein’s  dayside  aurora  and  this 
one  is  the  fact  that  the  first  one  covered  only  winter  arcs  at  night,  whereas  our 
work  includes  daylight  aurora  that  coexists  with  dayglow.  Although  background 
dayglow  has  been  subtracted,  the  physics  of  the  emitting  species  (e.g.  density, 
height,  etc.  )  and  atmospheric  transfer  function  may  differ. 

Also  consider  the  situation  of  a  temporally  dynamic  dayside  aurora,  even  at 
low  tfps.  Averaging  two  small  samples  of  arcs  may  result  in  large  differences  in 
average  location  and  intensity. 
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Figure  .21:  Auroral  frequency  -  Kp  =  2,  A  —  1 356 A . 


42 


Figure  .22:  Auroral  frequency  -  Kp  =  3,  A  =  1596A. 
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Figure  .24:  Auroral  frequency  -  Kp  =  4,  A  =  1596-4. 
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Figure  .29:  Average  auroral  intensity  -  Kp  =  1,  A  =  1356.A. 
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Figure  .33:  Average  auroral  intensity  -  Kp  —  3,  A  =  1356^. 
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.7.2  Future  Research 


This  report  summarises  the  first  two  years  of  research  on  image  processing  of 
the  Polar  BEAR  data.  The  software  has  matured  to  a  complete  package  for  the 
pre-processing,  geometric  and  photometric  calibration,  and  post-processing  of 
AIRS  images. 

Figures  16-35  of  this  report  depict  first  results  of  a  project  recently  initiated 
utilizing  the  above  package.  Our  aim  is  to  create  a  complete  and  comprehensive 
set  of  maps  of  auroral  frequency  and  average  brightness  in  the  far  ultraviolet 
as  a  function  of  geomagnetic  activity  indices.  The  semiautomatic  nature  of 
the  image-processing  routines  will  allow  us  to  enlarge  the  small  sample  (~  40 
scenes)  to  create  a  statistically  valid  database  (~  50  scenes  per  Kp  group)  for 
the  above  maps. 

We  will  correlate  known  geomagnetic  activity  indices,  such  as  IMF  and 
Kp,  with  the  well-resolved  discrete  arcs  which  are  immersed  in  the  continuous 
aurora.  The  sensitivity  of  at  least  one  of  the  detectors  (1356A)  affords  us  to 
study  and  compare  the  magnetic  latitudinal  and  longitudinal  distribution  of 
the  localized  arcs  with  the  continuous  aurora.  Basic  attributes  of  the  aurora, 
such  as  intensity,  dimensional  morphology,  and  occurrence  probability,  will  be 
quantitatively  examined. 

Recently,  the  daylight  aurora  has  captured  the  imagination  of  many  re¬ 
searchers.  Polar  BEAR  images,  processed  with  the  package  we  have  developed, 
present  an  unmatchable  opportunity  for  the  study  of  this  unique  and  unexplored 
territory.  The  airglow  subtraction  technique  allows  one  to  study  daylight  aurora 
without  the  interfering  background  radiation.  Furthermore,  relations  between 
airglow  indices  such  as  solar  flux,  geomagnetic  coordinates,  and  electron  densi¬ 
ties  may  permit  true  geographic-coordinate  mapping  of  electron  density  profiles. 
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